New analytical progress in the theory of vesicles under linear flow 
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Vesicles are becoming a quite popular model for the study of red blood cells (RBCs). This is a 
free boundary problem which is rather difficult to handle theoretically. Quantitative computational 
approaches constitute also a challenge. In addition, with numerical studies, it is not easy to scan 
within a reasonable time the whole parameter space. Therefore, having quantitative analytical 
results is an essential advance that provides deeper understanding of observed features and can 
be used to accompany and possibly guide further numerical development. In this paper shape 
evolution equations for a vesicle in a shear flow are derived analytically with precision being cubic 
(which is quadratic in previous theories) with regard to the deformation of the vesicle relative to 
a spherical shape. The phase diagram distinguishing regions of parameters where different types 
of motion (tank-treading, tumbling and vacillating-breathing) are manifested is presented. This 
theory reveals unsuspected features: including higher order terms and harmonics (even if they 
are not directly excited by the shear flow) is necessary, whatever the shape is close to a sphere. 
Not only does this theory cure a quite large quantitative discrepancy between previous theories 
and recent experiments and numerical studies, but also it reveals a new phenomenon: the VB 
mode band in parameter space, which is believed to saturate after a moderate shear rate, exhibits 
a striking widening beyond a critical shear rate. The widening results from excitation of fourth 
order harmonic. The obtained phase diagram is in a remarkably good agreement with recent three 
dimensional numerical simulations based on the boundary integral formulation. Comparison of our 
results with experiments is systematically made. 

PACS numbers: 87.16.D- 83.50.Ha 87.17.Jj 83.80.Lz 87.19.rh 



I. INTRODUCTION 

A vesicle is a closed simply connected membrane sepa- 
rating two liquids. Its dynamics in an ambient flow is an 
important research target due to various medical appli- 
cations. For instance, they are employed as biochemical 
reactors 0| , as vectors for targeted drug and gene de- 
livery [1,1J|, and as artificial cells for hemoglobin encap- 
sulation and oxygen transport Q. 

One can also regard a vesicle, from mechanical point 
of view, as a simplified model of a living cell. The most 
prominent biological counterpart for a vesicle is the red 
blood cell (RBC). The vesicle system allows to iden- 
tify the elementary processes of visco-elastic properties 
of cells moving passively in a fluid. The study of a single 
vesicle under a shear flow is one of the simplest nonequi- 
librium example. Nevertheless, this problem gives rise to 
a very complicated dynamics with different types of vesi- 
cle motion because of the nontrivial competition between 
the interactions of straining and rotational parts of shear 
flow with the vesicle. Three main types of motion for a 
vesicle in a shear flow have been identified so far: (i) tank- 
treading (TT), when shape and orientation of the vesicle 
is a steady-state, (ii) tumbling (TB), when the vesicle 
makes full rotations, and (iii) vacillating-breathing (VB, 
aka trembling or swinging) , when the longest axis of the 
vesicle oscillates around certain direction with the shape 
strongly changing during these oscillation cycles. Since 
many properties of a vesicle (e.g. the effective viscosity) 



strongly depend on its type of motion, it is important 
to define the phase diagram of a vesicle in a shear flow, 
predicting which type of motion will be realized for each 
set of parameters defining the vesicle dynamics. 

A widely accepted model for a vesicle proposes the fol- 
lowing assumptions: zero Reynolds number limit, mem- 
brane impermeability and local inextensibility, and the 
continuity of velocity field and mechanical stress across 
the membrane. The force exerted on the membrane by 
the liquids is balanced by the membrane rigidity force 
that is calculated from the energy of membrane bending 
taken aslol 



= y" {2H fdA + J ZdA, 



(1) 



where n is the bending rigidity coefficient of the mem- 
brane, H is the mean curvature and Z is a Lagrange 
multiplier, ensuring the local membrane inextensibility. 
Under these assumptions, and after a certain choice of 
units for all variables, the vesicle dynamics can be com- 
pletely described by three dimensionless numbers (see for 
example Q): the viscosity contrast 
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the excess area relative to a sphere 



A = Ar^'^ - 47r, 



(2) 



(3) 
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and the capillary number 



Ca 



(4) 



Here r]int,ext are the viscosities of fluids inside and outside 
the membrane respectively, A is the surface area of the 

1/3 

vesicle, = (f^) is the radius of a sphere containing 
the same volume as the vesicle, and 7 is the shear rate. 
It is convenient to choose ro as the unit of distance so 
that the volume of the vesicle is 
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(5) 



We also use T]ext as the unit of viscosity and 7^^ as the 
unit of time. 

Several studies have been recently dedicated to the 
determination of the phase diagram regarding different 
types of motion of an almost spherical vesicle in a shear 
flow. These studies were theoretical (8l-fll| . experimental 
[]"3 - tl5| , and numerical [TgI . [17| . Those theoretical works 
treat the membrane deviation from a spherical shape as 
a perturbation (e = -v/A is the small expansion parame- 
ter), represented as a series of spherical harmonics. Then 
it is suggested that one or two first orders in the expan- 
sion of the shape evolution equations must be retained. 
Generally, the resulting equations differ only because the 
authors propose distinct rules of neglection. The Leading 
Order TheoryQ (LOT) keeps only the terms of order of 
the imposed flow and the effects of the vorticity in the 
final equations. This theory is precise up to 0{e) in the 
expansion scheme. Later studies have included higher or- 
der terms in the expansion (up to 0(e^)). Lebedev et al. 
[Tl| included the next order terms in the membrane Hel- 
frich force (as compared to LOT), while another study 
(hereafter called Higher Order Theory[l3| (HOT)) ac- 
counted for the next order terms not only in the Helfrich 
force but also in the hydrodynamic field as well. 

In the theories of Ref . and Ref . [ll| , as well as in 
numerical studies based on dissipative particles dynamics 
[l7| . a saturation of the VB/TB phase border for large 
enough Ca was suggested- i.e. the value of A at which the 
transition from VB to tumbling occurs becomes almost 
independent on the capillary number (or equivalently on 
shear rate). 

Recent advances in three dimensional numerical com- 
putations based on the boundary integral formulation 
have made it possible to study vesicle dynamics 
quantitatively. The results of this study show, that the 
loss of stability of the TT solution occurs at values of 
A significantly higher than those predicted by analytical 
calculations [13, Furthermore, no saturation of the 
VB/TB phase border upon increasing shear rate (or Ca) 
was found. These results provide a new input that is 
worth understanding from the analytical point view. 

The main objective of this work is to investigate the 
reasons for the discrepancies between the results provided 
by the recent numerical simulations jT^ and analytical 



theories. It was pointed out in Ref.[ll| that the critical 
values of A at which the steady-state solution loses its 
stability, scales as 0(A~^/^) for a fixed shear rate. How- 
ever, this behaviour, which also agrees with the result 
of Keller and Skalak [l^, does not provide the correct 
information about the next order terms, as will be seen 
here. The next correction turns out to cause a rather big 
shift of the phase borders. 

This peculiar behavior (i.e. that the critical A diverges 
at vanishing A) confers to the vesicle problem a special 
status: in the small A limit, one further order and only 
one in the expansion scheme in powers of A is needed 
as compared to previous studies. As a consequence, it 
will be shown, in particular, that the next order term, 
previously unaccounted for, survives whatever small the 
deviation from a sphere is. This is the source of devia- 
tion between the recent numerical work and the existing 
analytical theories. As will be seen inclusion of one more 
order in the expansion will allow us to extract a phase 
diagram that is in quantitative agreement with the re- 
sults of numerical simulations [l6| . It will also be shown 
that fourth order spherical harmonics (in previous theo- 
ries only second order harmonics were included) can not 
be neglected, and especially beyond a certain shear rate. 
Indeed, it will be revealed that these harmonics give rise 
to a new important feature recently revealed in numerical 
simulations [l^. More precisely, the VB/TB phase bor- 
der does not saturate for fixed A when shear rate exceeds 
a certain value, rather a significant widening is revealed, 
in contrast to previous theoretica l fiol . [ll| , numerical 
and experimental investigations [l4| . 



II. SHAPE EVOLUTION EQUATIONS 
A. Basic definitions and the expansion scheme 

Here we perform the expansion of the shape evolution 
equations one order higher as compared to the most ad- 
vanced of previous papers [T^- Since a slightly deflated 
vesicle in the absence of flow takes a centro-symmetric 
shape, and the shear flow has a symmetry center in ev- 
ery point, we consider only symmetric shapes. 

Naturally, we take the origin of coordinates at the cen- 
ter of the vesicle moving with the velocity of the undis- 
turbed shear flow. That the vesicle moves with the same 
velocity as the applied shear flow is a consequence of 
the symmetry of the problem. We use the conventional 
parametrization of the vesicle membrane with reference 
points belonging to a sphere of radius unity 



R = x[l + f{x)]. 



(6) 



The shape function f{x) is then expanded in spherical 
harmonics of x. Note that the shape depends also on 
time, but the temporal variable will be specifled only 
when need be. Following Ref.jl, [l^ we introduce a for- 
mal expansion parameter e = \/A which helps classifying 
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different orders. Only second order harmonics are present 
to the order of 0{e) in the equilibrium shape function of 
an almost spherical vesicle in the absence of flow. They 
are also the only harmonics that can interact with linear 
flow directly. For consistency considerations, pushing the 
expansion to next order implies that the fourth and ze- 
roth order harmonics cannot be neglected. They enter 
dynamics as a result of interactions between the vesicle 
shape and the Helfrich force and the flow; they are of or- 
der 0{e^). The shape function does not contain spherical 
harmonics of odd orders owing to the centro-symmetry. 
The spherical harmonics of even orders higher than 4 
have an amplitude 0{£^) and thus are neglected. As fol- 
lows from the inequality 



1 = ^^ 
An 



An 



the average radius of the vesicle Rq is not exactly 1. 
Therefore a negative correction must be added to it in 
order to satisfy ([5]). This correction is of order O(e^) 
and leads to a non-zero coefficient for zero-order spher- 
ical harmonic in the expansion of f{x), which will be 
denoted /o below. 

The expansion takes formally the form: 

f{x)^e\fo + ef2 + e^f4 + 0{e^), (7) 



l=-k 

There are various definitions of spherical harmonics, 
here we use those that satisfy the following normalization 
conditions: 



Yuj{^)Yli{x)d?x 



47r(fc!) 



{2k + l){k + l)\{k-l)\ 



However, all equations are written in universal form and 

are valid for any rescaling of spherical harmonics. 

The condition together with the definition A = 
provide two additional constraints on the coefficients 

of the expansion We use the former to express /o 

through /2,i 

foix)d'x^-J h{xfd^x-'- J h{xf<fx + 0{e^) 

(8) 

and the latter to norm the coefficients f2^i 

2 j f2{xfdPx - y 1 UxfdPx = 1 + 0{e^). (9) 

We retain one more order in these expressions compared 
to previous works for the sake of consistency. Note that 
the volume and the surface have been expanded to the 
order O(e^) in order to obtain the equalities ([5]) and 



B. Derivation of the evolution equations of the fh 
from the boundary integral formulation 

We shall adopt here a different spirit from that of pre- 
vious theories |8l-[ll|. There the Lamb solution is used 
(a solution of the Stokes equations inside and outside 
the vesicle). We find it convenient to take the bound- 
ary integral formulation as a starting point. The present 
spirit is equivalent to using the Lamb solution, but it has 
the advantage that the boundary conditions at the mem- 
brane are already implemented in the boundary integral 
formulation. Indeed, the Stokes equations together with 
the boundary conditions on the membrane, and far away 
from it, can be converted into boundary integral formu- 
lation [131 . This leads to the following integral equation 

vi{x){l + \) = 

= 2ui{x) + 2 / Gji {R{x) - R{x')) Fj{x')d^R{x')+ 



(1 



A) J Kjira {R{x) - R{x')) vAx')Nm{x')d^R{x'). 



(10) 



Here v(x) is the actual velocity of point R{x), u{x) is 
the velocity of point R{x) in the imposed fiow, F{x) 
is the membrane rigidity force at point R{x), N{x) is 
the outward pointing normal to the vesicle surface at 
point R{x). The integration is taken over the surface of 
the vesicle and d^R{x') is the surface area element. The 
integration kernels have the following form: 



Gij {R) 
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Using the expression for the Helfrich force, we can find 
from (fTO|) the velocity field on the surface of the vesicle. 
The Helfrich force is given in terms of the shape as 



F^=-2 {k[2H{H^ - K 
dZ{x) 



AsH] - ZH) Nr 



dRj{x) 



(11) 



where K is the Gaussian curvature and A5 is the 
Laplace-Beltrami operator. It can easily be checked that 
expression (jlip vanishes for a spherical shape, and thus 
is of order 0(e). In order to balance pil|) at order 0(1) 
we need to assume the imposed flow to be of the same 
order as (jlip . This requirement can equivalently be ful- 
fllled, following Refs. H, [l3|, by the demand that n scale 
as e^^. We then set k = with R of order 0(1). 

The precise technical details will be given elsewhere, 
while here we shall present only the spirit and some in- 
termediate steps. We expand v{x), u{x), and F{x) into 
vector spherical harmonics of x. For convenience we de- 
flne them as 



Ylk.iix) 



{2k + l)x^Yk,l{x) ~ d,Yk,iix), 
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3,k.l\ 



Here the differentiation with respect to Xi is taken for- 
mally as if X were a regular 3D vector not bound to a 
sphere of radius unity, and eimn is the Levi-Civita (unit 
anti-symmetric) tensor. The advantage of such a defini- 
tion is that not only do these spherical harmonics consti- 
tute an orthogonal set with respect to the integration of 
dot product over x, i.e. the quantity 



is zero unless j = j' , k = k', I = I', but also operators 
G and K are diagonal in the chosen basis for a spherical 
vesicle, so that the integrals 



Ylk.iix)Gim{x — x')Y"[i^, J, (x'y xd^ x\ 



Yj,k.l{-^)-l^imn{x X )YJl^j^, J, [x ) Xnd'^xd'^X 



are zero unless j ~ j', k ^ k', I = V . Given the centro- 
symmetry we impose, in our expansion of t;(a;), u{x)^ and 
F{x), the absence of Y^ j, i{x), for even fc, and Y^ i{x), 
and Yr^ f, i{x), for odd k. In order to find the evolution 
equations we need to expand the velocity field to the 
order of O(e^). Coefficients for Y^^^iix) for fc > 6 are 
0{e^) and thus can be neglected. It will be shown later 
that coefficients of 5th and 6th orders of vector spherical 
harmonics make corrections to the evolution equations 
for fi^ i with fc < 4 which are O(e^). So only five first 
orders (starting from the zeroth one) of vector spherical 
harmonics should be taken into account in the expansion 
of any space-dependent quantity under consideration. 

In order to fulfil the condition of local inextensibility 
of the membrane we impose zero surface divergence of 
the velocity field Q. 



dA{x) dvi{x) 



dt 



dRj{x) 



[% - N,{x)Nj{x)] dA{x) = 0. (12) 



Here 



dvi{x) dvi{x) dxk 



dRj{x) dxk dRj{x) 

is the Jacobian matrix. We can take the derivatives as 
if X were a regular 3D vector, because the surface di- 
vergence is fully determined by the distribution of the 
velocity field on the membrane and does not depend on 
the continuation of the flow into the liquids. 

We expand equations (|10ll2p to order O(e^) and the 
resulting integrals could be performed analytically. Pro- 
jecting the results of the integration on the space of vector 
spherical harmonics up to the fourth order, we obtain a 
set of equations satisfled by the coefficients entering the 
expansions of v{x) and Z{x). Since the surface area is 



evaluated up to order O(e^), while the velocity field in 
(fT2| is expanded only up to order 0(£^), it is not appro- 
priate to use (fT2|) in order to ensure the conservation of 
the whole surface area. Therefore, we only use the projec- 
tions of on Yk.i{x) for k equal to 2 or 4. Accordingly, 
we leave the isotropic part of Z{x) (denoted here as Zq) 
undetermined. Once the final shape evolution equations 
are obtained we shall use the constraint that the time 
derivative of ^ is equal to zero in order to determine 
Zn. This way of reasoning was used in previous studies 
(sl-ITH. Note that unlike other theories [l-[ll|, since we 
expand the equations to higher order, it is not legitimate 
to replace dA{x) in p2|) with d'^x prior to projection 
on the spherical harmonics sub-space. Such a substitu- 
tion would imply neglecting terms of order O(e^) in final 
equations, and would be inconsistent with the spirit of 
the present theory. 

Having determined the velocity field thanks to the 
above expansions, we are in a position to obtain the fi- 
nal evolution equation by making use of the kinematic 
equation expressing the fact that the membrane velocity 
is equal to the fiuid velocity at the membrane 



df{x) dif{x)vj{x){6ij ~ XiXj) 
= v,{x)xi 



dt 



Rix) 



(13) 



Then the task is to substitute the expanded velocity 
field into this equation (in terms of coefficients of the 
velocity field) and project the resulting expression onto 
the space of spherical harmonics of interest. This then 
leads to the determination of the evolution equations that 
must be satisfied by the shape coefficients fi^i. It can be 
shown that the following identities hold l^g i{x)xi = 0, 
Yi^e^i{x)xi = 7Ye,i{x), and Y^q i{x)x^ = 6Ye^i{x). There- 
fore the projection of the first part of the right hand 
side of ([T^ on the subspace of spherical harmonics of 
orders up to four does not depend on the coefficients of 
the vector spherical harmonics of orders five and six in 
the velocity field. Regarding the the second part, dif{x) 
is of order 0(e) and the coefficients for fifth and sixth 
orders of vector spherical harmonics in the velocity field 
expansion are of order 0{£^) so that their contribution 
in the shape evolution equation is of order O(e^), that is 
beyond our accuracy. 

It is convenient to decompose the applied shear flow 
into its elongational and rotational parts to simplify 
the flnal equations. The same decomposition can be 
applied to general linear flow: the quantity E2{x) = 
Ui{x)xi/2 defines the straining part, while the vector 
= eijkdjUk{x) /2 represents the vorticity. Note, that 
u{x) in pni) is the velocity of the imposed fiow evaluated 
at the point R{x), which takes the following form 

ui{x) = {diE2 + eijkXjnk)il + fix)) 

The shape evolution equations for the second and 
fourth order harmonics can be written in a compact form: 



Df2j ^ J F2{x)Y2.i{x)*d^x 
Dt jY2,i{x)Y2.i{x)*d^x' 



(14) 
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Dt jY4,i{x)Yij{x)*d^x' 



(15) 



The left hand side term of the equation is a special 
derivative of fk,i that naturally arises in a non-rotating 
coordinate system, and its definition is 



Dfk,i dfkj , / eij„id.ifkXjUJmYkj{x)*d'^x 



Dt 



dt 



jYk,i{x)Yk,i{x)*d^x 



(16) 



To the order at which our expansion is performed, the 
rotational velocity of the vesicle is not equal to the vor- 
ticity of the imposed flow (unlike lower order calculations 
[SHllI), but rather it has a new contribution originating 
from the shape function 



ejkidikE2diif2 



(17) 



Note that if only the first term Hi is retained (8l-[ll| , then 
(|16p coincides with the Jaumann derivative. 
The functions Fk can be written as 

F2 = aif2+ea2fi+e\a3f^ + a4f2f4) + biE2+eb2E2f2+ 
+e\b3f2d,f2d^E2 + bidj2d,,f2d,E2 + 
+b5d^f2d^jE2djf2 + bQE2U) + O(e'), 

F4 = £(ci/4 + C2/I) + £2(C3/| + C4/2/4) + ediE2f2 + 

+e'{d2hE2 + d3f^E2 + d5f2d,E2d,f2) + Oie^). 

The coefficients a^, bi, Ci^ and di are rational functions of 
A, a = e/Ca, and Zq. The exact expressions are listed in 
the appendix. 



III. THE PHASE DIAGRAM 

A. The analytical phase diagram and comparison 
with previous works 

The phase diagram for A = 0.43 is presented on the 
Fig.l. The value A = 0.43 is chosen for the sake of com- 
parison with available numerical data [l6| . Also shown 
on that Fig ure are the results obtained in two previous 
theories [id, [U. As can be seen the basic structure of 
the phase diagram bears similarity with previous ones. 
There are however some important distinctions. First, 
the phase borders are significantly shifted towards the 
region of higher viscosity contrasts. Second, the VB/TB 
transition curve does not saturate upon increasing shear 
rate (or Ca), rather it exhibits a striking widening. We 
have also compared the results with those obtained re- 
cently by full three dimensional simulations. A remark- 
able agreement between the two approaches is found (see 
comparison in Ref. (l6|). 




Figure 1: Continuous lines represent the phase diagram for 
A — 0.43, the results of the higher order theory and Lebe- 
dev et al. theory [iJl are added in dashed and dotted lines for 
comparison. The TT phase borders are almost indistinguish- 
able for last two theories. 



B. Basic reasons for necessity of higher order 
expansion 

Our first concern is to understand why the inclusion of 
the terms of order 0{e^) provides such a dramatic shift of 
the phase borders even though the shape is almost spher- 
ical (the relative excess area is only 0.43/(47r) = 0.034). 
We have studied the evolution of the phase diagram as a 
function of A in order to gain further insight. We have 
tracked the viscosity contrast A at which the loss of sta- 
bility of TT motion occurs in the two asymptotic limits 
— > (the corresponding value is denoted as Acq) and 
Ca ~ CO (denoted as Acoo)- It was reported in [11| that 
Ac oc . Since we expect the expansions of Ac to be an- 
alytic in e, we have thus attempted the following ansatz: 



Ac = A(-i)e-i + Ai°) + 0(£). 



(18) 



To check the validity of this expansion we plot these crit- 
ical values as a function of A^^/'^ = 1/e. We expect then 
a linear behavior. This is presented on Figj^] where we 
also compare our results with those of previous theories. 
A key point to be discussed further below is the fact that 
previous theories provide correct values for the dominant 
term Ac in the expansions (|18p . but not for A^^' which 



does not vanish even in the limit of almost spherical vesi- 
cles (note that the next order terms in the expansion pS)) 

tend to with the excess area). Because A^'^'' remains fi- 
nite whatever small is the deviation from a sphere, the 
discrepancy between the TT phase borders obtained by 
different theories (with one lower order below the present 
one, namely of the order of 0{e)) persists for any excess 
area. In contrast the present theory valid to order O(e^), 
has the property that even if one wishes to make an ex- 
pansion to the next order (i.e. order 0(e^)), then the 
shift in the borders in the phase diagram would be negli- 



6 




Figure 2: Ao and Xoo as a function of A 



gibly small, i.e. the shift would vanish in the small excess 
area limit. In other words, the present theory shows that 
there is a convergence of the small deformation scheme at 
order O(e^). This implies that there is no need to continue 
the expansion of the shape evolution equation beyond the 
order O(e^) since this would neither notably modify the 
results for A < 0.43 nor significantly extend the applica- 
bility to higher values of A. Comparison of the present 
theory with numerical results (see comparison provided 
in Ref. [l^) shows a good agreement for A = 0.43. (the 
only value explored to date by the numerical scheme). 

C. Discussion of the tank-treading phase borders 

As can be seen the dependencies shown on Fig|2] are 
fairly linear for A < 1 and are almost indistinguishable 
from straight lines for A < 0.5. The region of l/-\/A 
where the curves noticeably deviate from straight lines 
is wider for present theory than for HOT. This is most 
likely caused by the introduction of fourth order spherical 
harmonics, their amplitude is no more small compared to 
that of the second order harmonics for A > 1. The ex- 
pression ([9]) becomes inapplicable, and as a result fi,i{t) 
grow uncontrollably with time in solutions of the differ- 
ential equations. This problem can be fixed by expanding 
the excess area to the order 0{e^), but the results of such 
patching are not trustworthy because they go beyond the 
precision of the differential equations. 

Note also that the second term in (fT7|) tends to align 
the vesicle to the elongating direction of the straining 
part of the shear flow which makes the angle 7r/4 with 
the direction of the flow. In the TT phase close to the 
loss of the stability of the steady-state solution the vesi- 
cle is almost parallel to the flow. The effective vorticity 
([17]) turns out to be less than in the TT phase, and 
it leads to the increased stability of the steady-state so- 
lution. The correction to the vorticity is proportional to 
£ and thus increases with e and dominates over in the 



Theory 






LOT[8] 


|;V307i^(3.38) 


-§(-1.39) 


HOT[10| 


3.91 


-1.22 


Lebedev et a/.[lll 


if VT07r(3.90) 


-§(-1.39) 


Keller and Skalak [181 


2\/TC)^(3.74) 


1|(1.17) 


Present 


3.89 


1.54 



Table I: X).^ and AJ^q extracted from different theories. It 
should be noted that Lebedev et al.[lll| theory provided only 
values for A^g^' and we have formally extracted from their 
theory the value of A[,g' for comparison purpose. 



Theory 




i(0) 


LOT[8j 


|;V30^(3.38) 


-§(-1.39) 


HOT[10J 


4.79 


-1.38 


Lebedev et a/.[lll 


i| 715^(4.78) 


-§(-1.39) 


Present 


4.76 


1.13 



Table II: Aico ' and aISj extracted from different theories. It 
should be noted that Lebedev et al.[lH theory provided only 
values for A^^^' and we have formally extracted from their 
theory the value of X^ciL for comparison purpose. 



large excess area limit. As a consequence, the truncation 
PT)) becomes illegitimate. Generally, we may assert that 
A < 1 is a good estimate for the applicability limits of 
the small deformation approximation. 

Some theories[l, [ll|, [l^ allow for analytical extraction 
of the coefficients Ai and x'^cK and for HOT and the 
present theory the extraction is numerical by using the 
slope and offset of the lines on the Fig. 2. The combined 
data are presented in the tables 1 and 2. 

HOT and Lebedev et al. theories show a good agree- 
ment for the coefficients Ac"^' with the present calcula- 
tion. This means, that expansion of the shape evolution 
equations to the order 0(e) was indeed sufficient in or- 
der to capture the correct value of even for large 
values of Ca, where Lebedev et al. theory loses its ap- 
plicability. The present theory shows that the correct 
value of X^c^ is captured thanks to the fact that the ex- 
pansion scheme is pushed one step further in the excess 
area. Keller-Skalak theory gives a slightly different result 
for aI , the difference is caused by the fact they used a 
velocity field having a surface divergence of order 0(e), 
whereas the membrane local inextensibility requires zero 
surface divergence instead. On the other hand Keller- 
Skalak theory provides a better estimate for the value of 
X^c^ than the one proposed by HOT and this explains why 
its results are relatively close to the results of numerical 
simulations for A ~ 1. 
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D. Discussion of TB/VB transition 



We plan to dedicate a separate research to the proper- 
ties of different types of vesicle motion, so here we only 
briefly discuss further implication of the new theory on 
the TB/VB transition. Neglecting harmonics of orders 
higher than two, as done in previous theories, provides 
two major simplifications. First, the vesicle shape has 
three symmetry planes, which make the definition of the 
orientation angle obvious, and second the in-plane mo- 
tion is defined by only two independent variables. As a 
consequence, the dynamics relaxes with time to a cyclic 
motion (which degenerates into a point in the TT phase). 
Thus only a simple limit cycle can exist, and this shows 
that each of the VB and TB modes has its own region of 
existence in parameter space. The above assertions can 
not be made in general, and especially when the fourth 
order harmonics are included. It also turns out that near 
the transition region the vesicle tends to finish its TB 
quasi-cycles by assuming an oblate almost axisymmetric 
shape in the shear plane rather than an elongated shape 
perpendicular to the flow. For such shape the error in 
the definition of the orientation angle is very large. 

Let us now discuss how each type of motion is de- 
termined from our evolution equations. We start with 
some initial values for fki{t) and wait a ceratin time in- 
terval until the initial data are irrelevant. Then it was 
checked whether during one oscillation quasi-cycle f{x) 
has a maximum for x lying in the shear plane and perpen- 
dicular to the shear velocity. Despite some noise due to 
the aforementioned complications it is clearly seen that 
unlike with previous theories the VB/TB phase border 
does not saturate even for Ca = 10, and the VB region 
broadens with the increase of the capillary number (Fig. 

Note that it was suggested recently that higher 
orders of spherical harmonics can be excited close to the 
VB/TB phase border and this may cause some widen- 
ing of the VB phase region (but still saturation at large 
Ca is found, unlike our theory). It can be checked that 
when Zq + 20e/Ca < 0, the fourth order harmonics are 
excited (i.e. the corresponding decay time becomes infi- 
nite to leading order). Here, we found that close to the 
VB/TB transition at some times during the oscillation 
the quantity 



ZoCa 



(19) 



exceeds 20. This is probably the reason why the VB/TB 
phase border remains unsaturated for much larger values 
than in previous analytical studies. It was proposed in 
[20j . however, that excitation of higher harmonics was 
due to thermal fluctuations. This contrasts in spirit with 
our theory where the fourth order harmonic is excited as 
non-linear interaction of second harmonics and no refer- 
ence to temperature is required. 



A=0.16 
A=0.43 
A=0,8 

Lebedev er a!. Theory 




2 4 6 8 10 12 J4 16 18 20 



Figure 3: Phase diagrams for various values of A in 5 — A 
coordinates. 



IV. COMPARISON WITH EXPERIMENTS 

Experiments have been performed recently regarding 
determination of the phase diagram of vesicle motions 
under shear flow 



14| . They referred to a previous the- 



oretical work [ll| which had suggested that the phase 
diagram should depend on only two independent dimen- 
sionless control parameters (and not three), namely 



A 



i + |a|VA,5 



7^ Ca 



(20) 



3V3 A 

It has been reported [l3| (with a certain degree of uncer- 
tainty) that the experimental data were consistent with 
the fact that only the above two parameters determine 
the phase diagram. Our results do not comply with this 
report, as shown on Figure [3] Indeed, besides 5* and A, 
the excess area A plays an important role. Experiments 
mixed data for different A's in the plane (S*, A). The 
band of the VB mode in experiments looks quite wide, 
and we believe that this reflects the sensitivity of the lo- 
cation of the VB band to excess area (in other words 
the experimental data may be viewed as juxtaposition of 
bands each representing a value of A; see the example 
of Figure [31) . It is hoped that a study representing each 
value of A will be performed in the future with the aim 
of making comparison with theory clearer. 

Furthermore, it must be noted that our estimates for 
the TT-VB phase borders are higher than those observed 
in experiments. This is not very surprising since tran- 
sients are found to be very long close to TT-VB transi- 
tion (as also discussed in our recent full numerical simu- 
lations (3), and therefore a flrm conclusion on the type 
of motion can not be made on the basis of the current 
experimental data, which are limited to only few periods 
of oscillation (while transient can exhibit up to hundreds 
of cycles as discussed recently [l^). In other words, this 
kind of long relaxation would convey the impression that 
dynamics is of VB type, whereas in reality it is a TT one. 
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Present Theory, A=0.43 

Present Theory, A=l 

Lebedev ff (?/. TTieory 




Figure 4: Phase diagrams of vesicle motions under general 
planar linear flows. The viscosity contrast A is fixed to 1 and 
the ratio of rotating and straining parts of the flow ujo/sq is 
varied. The theory of Lebedev et a/.[ll| is plotted in dotted 
lines for comparison. 



A=0.43, s =1/2 



Present Theory, 0)^=1/2 (simple shear flow) 

Present Theory, (0^=3/4 

Present Theory, O)^ =5/4 

Lebedev ei al. Theory, arbitrary' (S3 




Figure 5: Phase diagrams of vesicles under general planar 
linear flows. The straining and rotational parts of the flow 
are denoted as sq and uio respectively. 



grams for various values of A is indeed not as striking as 
for the simple shear flow (Fig (3]): the TT phase borders 
fall close to each other and VB bands overlap for a wide 
range of A. At the same time, the VB/TB phase borders 
still vary significantly with the excess area. However, by 
fixing A and varying A we can produce series of phases 
diagrams in the S and A plane corresponding to different 
values of wq/so, (see FiglS]). Here we find the same kind 
of discrepancy as on Fig|3] . 

Finally, let us compare our results with other set of ex- 
periments in the TT regimefisj. We represent the vesicle 
orientation angle 0o under strong shear flow(Fig. 3) (for 
the sake of comparison with experiments which were per- 
formed at high shear rates). It can be checked that the 
exact form of the bending energy becomes insignificant 
under strong flows, because the dominant contribution to 
the membrane force comes form the tension part which 
enforces incompressibility of the membrane. Unlike for 
the VB/TB phase border, the dependence of (/>o on A 
quickly saturates with Ca- We take Ca = 100 in our cal- 
culation. We use the artificial parameter A ([^0]) instead 
of A in order to show the effect of O(e^) terms taken into 
account by the present calculation. Note that the results 
of Lebedev it et al. theory do not depend on the 
excess area A in this representation, unlike our results 
and those reported by experiments (see Fig. 5). The re- 
sults of HOT depend also on A, but that dependence is 
weak so that they hardly differ from theoretical results 
of Ref.[ll| (we do not plot them here in order to avoid 
encumbering of the figure). Experimental results (l3| are 
provided for four different values of A : 0.15, 0.24, 0.42, 
and 1.43. We exclude A = 1.43, since we do not expect 
our theory to be applicable to a such large value of excess 
area. As can be seen, theoretical results and experimen- 
tal ones show a rather good agreement provided that one 
is not close to the TT/VB border (i.e. if the angle is 
not too close to zero). The discrepancy for small incli- 
nation angles may be attributed to thermal fluctuations 
or interactions with walls. However, a systematic theo- 
retical study of these factors should be performed before 
drawing conclusive answers. 



Other experiments [15| were performed under linear 
flow with applied velocity field Ux = (so + ^o)y, Uy ~ 
(sq — uio)x, Uz = 0, which is different from a simple shear 
flow. The experiments were performed with no viscos- 
ity contrast and the ratio of the rotational (ujq) and the 
straining (sq) parts of the flow was varied instead. The 
results of the experiments were plotted in coordinates 

representing the generalization of ([20]). It was then 
claimed [15| that the phase diagram depends on these 
two parameters only (at least with no viscosity contrast). 
The Figd] simulates this experiment using the present 
theory. The discrepancy between resulting phase dia- 



V. DISCUSSION AND CONCLUSION 

Here we have presented a small deformation theory of a 
vesicle in shear flow keeping one more order in the expan- 
sion of the shape evolution equations than prior studies. 
We have conflrmed the leading term in the asymptotic 
expansion of critical viscosity ratios determined by previ- 
ous analytical studies, but we also were able to determine 
accurately the next term in the expansion. This term is 
constant and survives no matter how small the excess 
area is. Moreover, this theory provides a correction to 
the phase borders that is significant in a wide range of 
excess areas. Unlike previous analytical works, but in 
agreement with the results of the numerical simulations 
|16| . we observe an unsaturated growth of the VB/TB 
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TT inclination angle for C =100 



A Experiment, A=0.15 

□ Experiment, A=0.24 

9 Experiment, A=0.42 

— Present TTieory, A=0.15 

- ■ Preseni Theory, A=0.24 

Present Theory, A=0.42 

■ ■ ■ Lebedev el al. Theory 




Figure 6: inclination angle of vesicles in TT phase for large 

Ca. 



phase border in a quite large range of the capillary num- 
ber. 

We hope that this work will incite future experimental 
research. It will also be interesting to investigate ex- 
perimentally the presence of widening of the VB band 
as a function of shear rate. Finally, we did not in- 
clude third order spherical harmonics in the expansion 
of the shape function, restricting our calculations only 
to centro-symmetric vesicles. In the full numerical simu- 
lation [l6| . no symmetry restriction is imposed and so 
far no manifestation of the third order harmonic has 
been observed. This confers to our assumption of centro- 
symmetry a certain legitimacy. 
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Appendix A: The shape evolution equations and the 
expression of various coefficients in terms of physical 
parameters 

We remind that 

F2 = aif2+ea2fl+e\a3fi + aif2h) + hE2+eb2E2f2+ 

+e\b3f2dj2d,E2 + bidj2d,,f2djE2+ 



F4 = e(ci/4 + 02/2) + £'(C3/| + Cif2fi) + 



+ediE2f2 +e^{d2hE2+d3f^E2 + d5f2d^E2d,f2) + 0{e^). 
The coefficients are then written as 



ai = -24 



Zq + 6k 
23A + 32 



02 = 24 



(49A + 136)Zo + (432A -I- 1008)k 
(23A-f 32)2 



bi 



120 



23A + 32 



62 = 2400 



A-2 



(23A-^32)2 



-40 



241A + 344 
(23A-h32)2 



ci = -40 



Zo + 20R 
19A + 20 



C2 



16 (92 - 3A)Zo + (1822A + 3112)k 
y (19A-t-20)(23A-l-32) 



di 



20 1047A + 1072 
y (23A-l-32)(19A-h20) 



d2 = 10 



257A - 32 



(23A + 32)(19A-F20) 



d3 = - 
^ 6 



1 -125055A3 + 594716A2 + 1107168A + 392320 



(2A-f-5)(23A + 32)(19A + 20)2 



8 (2479595A3 + 6703156A2 + 18601472A-h 16622592^ 
°^~~T75 (23A + 32)3(19A + 20) °~ 

8 78181390A3 + 390845256A2 + 713624832A -I- 429094912 _ 
"175 (23A + 32)3(19A + 20) 
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(18335A2 + 57376A + 44224)Zo + (318896A2 + 1026064A + 809600)k 
(19A + 20)(23A + 32)2 

1 78420885A'^ + 815632786A3 + 21935721 12A2 + 1954954752A + 481607680 
~ 70 (2A + 5)(19A + 20)(23A + 32)3 

1 490295475A4 + 3878418742A3 + 9801602064A2 + 9403966464A + 2954997760 
630 (2A + 5)(19A + 20)(23A + 32)3 

1 136559745A'' + 385825454A3 + 141580368A2 - 26564352A + 98232320 
~ 315 (2A + 5)(19A + 20)(23A + 32)3 

4 2552583A3 + 1777744A2 - 913584A + 315392 
^^~~63 (23A + 32)2(19A + 20)2 °" 

16 32203871 A'V 1181 14286A2 + 154157080A + 69803008.^ 
~63 (23A + 32)2(19A + 20)2 

_ 2 (21383A2 + 63844A + 44928)Zo + 48(13411A2 + 35444A + 23040)k 
~ (23A + 32)(19A + 20)2 

1 41587815A4 + 95846332A3 - 88522016A2 - 295841536A - 152842240 
* " 42 (2A + 5)(23A + 32)2(19A + 20)2 
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